Cargar las librerías
library(tidyverse)
library(sf)
library(mapview)
library(GADMTools)
library(ggspatial)
library(leaflet)
library(leaflet.extras2)
library(spdep)
library(spatstat)
library(raster)
library(smacpod)
library(ggspatial)
library(DT)
Caso de estudio
Cargar los datos
covid <- read_csv(url("https://zenodo.org/record/4915889/files/covid19data.csv?download=1"))
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## FECHA_CORTE = col_date(format = ""),
## UUID = col_character(),
## DEPARTAMENTO = col_character(),
## PROVINCIA = col_character(),
## DISTRITO = col_character(),
## METODODX = col_character(),
## EDAD = col_double(),
## SEXO = col_character(),
## FECHA_RESULTADO = col_date(format = ""),
## rango_edad = col_character(),
## lon = col_double(),
## lat = col_double()
## )
covid_p <- covid %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326)
covid_p %>%
filter(FECHA_RESULTADO=="2020-12-11") %>%
ggplot() +
geom_sf()
m_p<-covid_p %>% filter(FECHA_RESULTADO == "2020-12-10") %>%
mapview(layer.name="puntos")
m_p
peru <- gadm_sf_loadCountries("PER", level=3)
lima_sf <- peru %>%
pluck("sf") %>%
# Filtramos los datos espaciales solo de Lima metropolitana
filter(NAME_2 == "Lima") %>%
# Editamos algunos errores en nuestros datos espaciales
mutate(NAME_3 = ifelse(NAME_3 == "Magdalena Vieja",
"Pueblo Libre", NAME_3))
mexico<-gadm_sf_loadCountries("MEX",level=2)
cdmx_sf<-mexico %>%
pluck("sf")
mexico
## $basename
## [1] "./"
##
## $sf
## Simple feature collection with 1854 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -118.3689 ymin: 14.53292 xmax: -86.71014 ymax: 32.71804
## Geodetic CRS: WGS 84
## First 10 features:
## ISO NAME_0 NAME_1 NAME_2 TYPE_2 ENGTYPE_2
## 3683 MEX Mexico Aguascalientes Aguascalientes Municipio Municipality
## 3665 MEX Mexico Aguascalientes Asientos Municipio Municipality
## 3687 MEX Mexico Aguascalientes Calvillo Municipio Municipality
## 3652 MEX Mexico Aguascalientes Cosío Municipio Municipality
## 3690 MEX Mexico Aguascalientes Jesús María Municipio Municipality
## 3675 MEX Mexico Aguascalientes Pabellón de Arteaga Municipio Municipality
## 3657 MEX Mexico Aguascalientes Rincón de Romos Municipio Municipality
## 3661 MEX Mexico Aguascalientes San José de Gracia Municipio Municipality
## 3663 MEX Mexico Aguascalientes Tepezalá Municipio Municipality
## 3336 MEX Mexico Baja California Ensenada Municipio Municipality
## geometry
## 3683 MULTIPOLYGON (((-102.448 21...
## 3665 MULTIPOLYGON (((-102.229 22...
## 3687 MULTIPOLYGON (((-102.613 21...
## 3652 MULTIPOLYGON (((-102.2537 2...
## 3690 MULTIPOLYGON (((-102.5781 2...
## 3675 MULTIPOLYGON (((-102.4122 2...
## 3657 MULTIPOLYGON (((-102.19 22....
## 3661 MULTIPOLYGON (((-102.4122 2...
## 3663 MULTIPOLYGON (((-102.2376 2...
## 3336 MULTIPOLYGON (((-114.7158 2...
##
## $level
## [1] 2
##
## $hasBGND
## [1] FALSE
##
## attr(,"class")
## [1] "gadm_sf"
covid_count_mx <- covid %>%
group_by(DISTRITO, FECHA_RESULTADO) %>%
summarise(casos = n()) %>%
ungroup() %>%
complete(FECHA_RESULTADO = seq.Date(min(FECHA_RESULTADO, na.rm =T), max(FECHA_RESULTADO, na.rm = T), by="day"),
nesting(DISTRITO), fill = list(n = 0))
## `summarise()` has grouped output by 'DISTRITO'. You can override using the `.groups` argument.
covid_count <- covid %>%
group_by(DISTRITO, FECHA_RESULTADO) %>%
summarise(casos = n()) %>%
ungroup() %>%
complete(FECHA_RESULTADO = seq.Date(min(FECHA_RESULTADO, na.rm =T),
max(FECHA_RESULTADO, na.rm = T),
by="day"),
nesting(DISTRITO), fill = list(n = 0))
## `summarise()` has grouped output by 'DISTRITO'. You can override using the `.groups` argument.
covid_sf <- lima_sf %>%
mutate(DISTRITO = toupper(NAME_3)) %>%
full_join(covid_count, by = "DISTRITO", "FECHA_RESULTADO")
class(covid_sf)
## [1] "sf" "data.frame"
covid_sf %>%
filter(FECHA_RESULTADO == "2020-12-11") %>%
ggplot() +
geom_sf()
m_sf <- covid_sf %>%
filter(FECHA_RESULTADO == "2020-12-11") %>%
mapview(layer.name = "distritos")
m_sf
Múltiples capas
ggplot() +
geom_sf(data = covid_sf %>%
filter(FECHA_RESULTADO == "2020-12-11")) +
geom_sf(data = covid_p %>%
filter(FECHA_RESULTADO == "2020-12-11"))
m_p +m_sf
Visualización de datos espaciales
Patrones de puntos
covid_p %>%
filter(FECHA_RESULTADO == "2020-12-11") %>%
ggplot() +
geom_sf(aes(col = SEXO), alpha = .2) +
facet_wrap(.~SEXO)
p1<-covid_p %>%
filter(FECHA_RESULTADO == "2020-12-11") %>%
mapview(layer.name = "points", zcol = "SEXO", burst = T)
p1
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:raster':
##
## select
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
Dos o más variables
covid_p %>%
filter(FECHA_RESULTADO == "2020-04-11" |
FECHA_RESULTADO == "2020-12-11") %>%
ggplot() +
geom_sf(aes(col = SEXO), alpha = .2) +
facet_grid(SEXO~FECHA_RESULTADO) +
guides(col = F)
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
m1 <- covid_p %>%
filter(FECHA_RESULTADO == "2020-04-11") %>%
mapview(zcol = "SEXO", layer.name = "2020-04-11 - SEXO")
m2 <- covid_p %>%
filter(FECHA_RESULTADO == "2020-12-11") %>%
mapview(zcol = "SEXO", layer.name = "2020-12-11 - SEXO")
m1 | m2
Composición
covid_p %>%
filter(FECHA_RESULTADO == "2020-12-11") %>%
ggplot() +
geom_sf(data = covid_sf) +
geom_sf(aes(col = EDAD), alpha = .2) +
scale_color_viridis_c(option = "B") +
annotation_scale() +
annotation_north_arrow(location = "tr",
style = north_arrow_nautical)+
theme_bw()
Datos en polígonos
Una variable
covid_sf %>%
filter(FECHA_RESULTADO == "2020-12-11") %>%
ggplot() +
geom_sf(aes(fill = casos))
covid_sf %>%
filter(FECHA_RESULTADO == "2020-12-11") %>%
mapview(layer.name = "casos", zcol = "casos")
Dos o más variable
covid_sf %>%
filter(FECHA_RESULTADO == "2020-04-11" |
FECHA_RESULTADO == "2020-12-11") %>%
ggplot() +
geom_sf(aes(fill = casos)) +
facet_grid(.~FECHA_RESULTADO)
d1 <- covid_sf %>%
filter(FECHA_RESULTADO == "2020-04-11") %>%
mapview(zcol = "casos", layer.name = "2020-04-11 - casos")
d2 <- covid_sf %>%
filter(FECHA_RESULTADO == "2020-12-11") %>%
mapview(zcol = "casos", layer.name = "2020-12-11 - casos")
d1 | d2
Composición
covid_sf %>%
filter(FECHA_RESULTADO == "2020-12-11") %>%
ggplot() +
geom_sf(aes(fill = casos)) +
scale_fill_viridis_c(option = "F", direction = -1) +
annotation_scale() +
annotation_north_arrow(location = "tr",
style = north_arrow_nautical)+
theme_void()
Variación espacial del riesgo
covid_subset <- covid %>%
filter(FECHA_RESULTADO == "2020-05-05")
covid_win <- owin(xrange = range(covid_subset$lon),
yrange = range(covid_subset$lat))
covid_ppp <- ppp(covid_subset$lon,
covid_subset$lat,
window = covid_win)
densidad_raster_cov <- raster(density(covid_ppp, bw.ppl),
crs = 4326) %>%
mask(lima_sf)
densidad_raster_cov %>%
mapview()
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
Detección de clústers
Datos puntuales
Estadísticas de escaneo espacial (Spatial Scan Statistics-SSS)
covid_subset_posi <- covid %>%
filter(FECHA_RESULTADO == "2020-05-05") %>%
mutate(positividad = ifelse(METODODX == "PCR", 1, 0))
covid_scan_ppp <- ppp(covid_subset_posi$lon,
covid_subset_posi$lat,
range(covid_subset_posi$lon),
range(covid_subset_posi$lat),
marks = as.factor(covid_subset_posi$positividad))
covid_scan_test <- spscan.test(covid_scan_ppp,
nsim = 49, case = 2,
maxd=.15, alpha = 0.05)
covid_scan_test
## $clusters
## $clusters[[1]]
## $clusters[[1]]$locids
## [1] 1103 1371 1092 1375 1381 1104 1380 1378 1091 1074 1097 1376 1089 1076 1373
## [16] 1379 1096 1087 1088 1095 2093 2122 2117 2102 2106 1105 2121 1107 605 2086
## [31] 1098 592 1101 2111 596 2096 1070 599 2120 597 971 2094 2119 607
##
## $clusters[[1]]$coords
## [,1] [,2]
## [1,] -77.0472 -12.11374
##
## $clusters[[1]]$r
## [1] 0.02709301
##
## $clusters[[1]]$pop
## [1] 44
##
## $clusters[[1]]$cases
## 110344
## 31
##
## $clusters[[1]]$expected
## [1] 11.17098
##
## $clusters[[1]]$smr
## 110344
## 2.775046
##
## $clusters[[1]]$rr
## 110344
## 2.873837
##
## $clusters[[1]]$propcases
## 110344
## 0.7045455
##
## $clusters[[1]]$loglikrat
## [1] 20.05829
##
## $clusters[[1]]$pvalue
## [1] 0.02
##
##
##
## $ppp
## Marked planar point pattern: 2316 points
## Multitype, with levels = 0, 1
## window: rectangle = [-77.18195, -76.64829] x [-12.468713, -11.634232] units
##
## attr(,"class")
## [1] "spscan"
# Construimos el centroide del clúster
cent <- tibble(lon = covid_scan_test$clusters[[1]]$coords[,1],
lat = covid_scan_test$clusters[[1]]$coords[,2]) %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326, remove = F)
# Construimos el área del clúster en base al radio
clust <- cent %>%
st_buffer(dist = covid_scan_test$clusters[[1]]$r)
cluster <- mapview(clust, alpha.regions = 0, color = "red")
points <- covid_subset_posi %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
mapview(zcol = "positividad", alpha.regions = .4, alpha = 0)
cluster + points
Datos agregados (en polígonos)
Autocprrelación espacial (global):Moran
covid_sf_subset <- covid_sf %>%
filter(FECHA_RESULTADO == "2020-05-05") %>%
mutate(casos = replace_na(casos, 0))
covid.nb <- poly2nb(covid_sf_subset, queen=TRUE,snap = 0.13)
covid.lw <- nb2listw(covid.nb, style="W", zero.policy=TRUE)
moran.test(covid_sf_subset$casos, covid.lw)
##
## Moran I test under randomisation
##
## data: covid_sf_subset$casos
## weights: covid.lw
##
## Moran I statistic standard deviate = 3.0548, p-value = 0.001126
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.090491242 -0.023809524 0.001399995
Autocorrelación espacial (local):Getis Ord
breaks <- c(-Inf, -1.96, 1.96, Inf)
labels <- c("Cold spot",
"Not significant",
"Hot spot")
covid_lg <- localG(covid_sf_subset$casos, covid.lw)
covid_sf_lisa<-covid_sf_subset %>%
mutate(cluster_lg=cut(covid_lg, include.lowest = TRUE,
breaks = breaks,
labels = labels))
covid_sf_lisa %>%
ggplot() +
geom_sf(aes(fill=cluster_lg)) +
scale_fill_brewer(name="Clúster",
palette = "RdBu", direction=-1) +
theme_bw()
